home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zgeesx.f < prev    next >
Text File  |  1996-07-19  |  13KB  |  372 lines
  1.       SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
  2.      $                   VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
  3.      $                   BWORK, INFO )
  4. *
  5. *  -- LAPACK driver routine (version 2.0) --
  6. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  7. *     Courant Institute, Argonne National Lab, and Rice University
  8. *     March 31, 1993
  9. *
  10. *     .. Scalar Arguments ..
  11.       CHARACTER          JOBVS, SENSE, SORT
  12.       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
  13.       DOUBLE PRECISION   RCONDE, RCONDV
  14. *     ..
  15. *     .. Array Arguments ..
  16.       LOGICAL            BWORK( * )
  17.       DOUBLE PRECISION   RWORK( * )
  18.       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
  19. *     ..
  20. *     .. Function Arguments ..
  21.       LOGICAL            SELECT
  22.       EXTERNAL           SELECT
  23. *     ..
  24. *
  25. *  Purpose
  26. *  =======
  27. *
  28. *  ZGEESX computes for an N-by-N complex nonsymmetric matrix A, the
  29. *  eigenvalues, the Schur form T, and, optionally, the matrix of Schur
  30. *  vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
  31. *
  32. *  Optionally, it also orders the eigenvalues on the diagonal of the
  33. *  Schur form so that selected eigenvalues are at the top left;
  34. *  computes a reciprocal condition number for the average of the
  35. *  selected eigenvalues (RCONDE); and computes a reciprocal condition
  36. *  number for the right invariant subspace corresponding to the
  37. *  selected eigenvalues (RCONDV).  The leading columns of Z form an
  38. *  orthonormal basis for this invariant subspace.
  39. *
  40. *  For further explanation of the reciprocal condition numbers RCONDE
  41. *  and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
  42. *  these quantities are called s and sep respectively).
  43. *
  44. *  A complex matrix is in Schur form if it is upper triangular.
  45. *
  46. *  Arguments
  47. *  =========
  48. *
  49. *  JOBVS   (input) CHARACTER*1
  50. *          = 'N': Schur vectors are not computed;
  51. *          = 'V': Schur vectors are computed.
  52. *
  53. *  SORT    (input) CHARACTER*1
  54. *          Specifies whether or not to order the eigenvalues on the
  55. *          diagonal of the Schur form.
  56. *          = 'N': Eigenvalues are not ordered;
  57. *          = 'S': Eigenvalues are ordered (see SELECT).
  58. *
  59. *  SELECT  (input) LOGICAL FUNCTION of one COMPLEX*16 argument
  60. *          SELECT must be declared EXTERNAL in the calling subroutine.
  61. *          If SORT = 'S', SELECT is used to select eigenvalues to order
  62. *          to the top left of the Schur form.
  63. *          If SORT = 'N', SELECT is not referenced.
  64. *          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
  65. *
  66. *  SENSE   (input) CHARACTER*1
  67. *          Determines which reciprocal condition numbers are computed.
  68. *          = 'N': None are computed;
  69. *          = 'E': Computed for average of selected eigenvalues only;
  70. *          = 'V': Computed for selected right invariant subspace only;
  71. *          = 'B': Computed for both.
  72. *          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
  73. *
  74. *  N       (input) INTEGER
  75. *          The order of the matrix A. N >= 0.
  76. *
  77. *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
  78. *          On entry, the N-by-N matrix A.
  79. *          On exit, A is overwritten by its Schur form T.
  80. *
  81. *  LDA     (input) INTEGER
  82. *          The leading dimension of the array A.  LDA >= max(1,N).
  83. *
  84. *  SDIM    (output) INTEGER
  85. *          If SORT = 'N', SDIM = 0.
  86. *          If SORT = 'S', SDIM = number of eigenvalues for which
  87. *                         SELECT is true.
  88. *
  89. *  W       (output) COMPLEX*16 array, dimension (N)
  90. *          W contains the computed eigenvalues, in the same order
  91. *          that they appear on the diagonal of the output Schur form T.
  92. *
  93. *  VS      (output) COMPLEX*16 array, dimension (LDVS,N)
  94. *          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
  95. *          vectors.
  96. *          If JOBVS = 'N', VS is not referenced.
  97. *
  98. *  LDVS    (input) INTEGER
  99. *          The leading dimension of the array VS.  LDVS >= 1, and if
  100. *          JOBVS = 'V', LDVS >= N.
  101. *
  102. *  RCONDE  (output) DOUBLE PRECISION
  103. *          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
  104. *          condition number for the average of the selected eigenvalues.
  105. *          Not referenced if SENSE = 'N' or 'V'.
  106. *
  107. *  RCONDV  (output) DOUBLE PRECISION
  108. *          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
  109. *          condition number for the selected right invariant subspace.
  110. *          Not referenced if SENSE = 'N' or 'E'.
  111. *
  112. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  113. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  114. *
  115. *  LWORK   (input) INTEGER
  116. *          The dimension of the array WORK.  LWORK >= max(1,2*N).
  117. *          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
  118. *          where SDIM is the number of selected eigenvalues computed by
  119. *          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2.
  120. *          For good performance, LWORK must generally be larger.
  121. *
  122. *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
  123. *
  124. *  BWORK   (workspace) LOGICAL array, dimension (N)
  125. *          Not referenced if SORT = 'N'.
  126. *
  127. *  INFO    (output) INTEGER
  128. *          = 0: successful exit
  129. *          < 0: if INFO = -i, the i-th argument had an illegal value.
  130. *          > 0: if INFO = i, and i is
  131. *             <= N: the QR algorithm failed to compute all the
  132. *                   eigenvalues; elements 1:ILO-1 and i+1:N of W
  133. *                   contain those eigenvalues which have converged; if
  134. *                   JOBVS = 'V', VS contains the transformation which
  135. *                   reduces A to its partially converged Schur form.
  136. *             = N+1: the eigenvalues could not be reordered because some
  137. *                   eigenvalues were too close to separate (the problem
  138. *                   is very ill-conditioned);
  139. *             = N+2: after reordering, roundoff changed values of some
  140. *                   complex eigenvalues so that leading eigenvalues in
  141. *                   the Schur form no longer satisfy SELECT=.TRUE.  This
  142. *                   could also be caused by underflow due to scaling.
  143. *
  144. *  =====================================================================
  145. *
  146. *     .. Parameters ..
  147.       DOUBLE PRECISION   ZERO, ONE
  148.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  149. *     ..
  150. *     .. Local Scalars ..
  151.       LOGICAL            SCALEA, WANTSB, WANTSE, WANTSN, WANTST, WANTSV,
  152.      $                   WANTVS
  153.       INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
  154.      $                   ITAU, IWRK, K, MAXB, MAXWRK, MINWRK
  155.       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SMLNUM
  156. *     ..
  157. *     .. Local Arrays ..
  158.       DOUBLE PRECISION   DUM( 1 )
  159. *     ..
  160. *     .. External Subroutines ..
  161.       EXTERNAL           DLABAD, DLASCL, XERBLA, ZCOPY, ZGEBAK, ZGEBAL,
  162.      $                   ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
  163. *     ..
  164. *     .. External Functions ..
  165.       LOGICAL            LSAME
  166.       INTEGER            ILAENV
  167.       DOUBLE PRECISION   DLAMCH, ZLANGE
  168.       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
  169. *     ..
  170. *     .. Intrinsic Functions ..
  171.       INTRINSIC          MAX, MIN, SQRT
  172. *     ..
  173. *     .. Executable Statements ..
  174. *
  175. *     Test the input arguments
  176. *
  177.       INFO = 0
  178.       WANTVS = LSAME( JOBVS, 'V' )
  179.       WANTST = LSAME( SORT, 'S' )
  180.       WANTSN = LSAME( SENSE, 'N' )
  181.       WANTSE = LSAME( SENSE, 'E' )
  182.       WANTSV = LSAME( SENSE, 'V' )
  183.       WANTSB = LSAME( SENSE, 'B' )
  184.       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
  185.          INFO = -1
  186.       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
  187.          INFO = -2
  188.       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
  189.      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
  190.          INFO = -4
  191.       ELSE IF( N.LT.0 ) THEN
  192.          INFO = -5
  193.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  194.          INFO = -7
  195.       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
  196.          INFO = -11
  197.       END IF
  198. *
  199. *     Compute workspace
  200. *      (Note: Comments in the code beginning "Workspace:" describe the
  201. *       minimal amount of real workspace needed at that point in the
  202. *       code, as well as the preferred amount for good performance.
  203. *       CWorkspace refers to complex workspace, and RWorkspace to real
  204. *       workspace. NB refers to the optimal block size for the
  205. *       immediately following subroutine, as returned by ILAENV.
  206. *       HSWORK refers to the workspace preferred by ZHSEQR, as
  207. *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
  208. *       the worst case.
  209. *       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
  210. *       depends on SDIM, which is computed by the routine ZTRSEN later
  211. *       in the code.)
  212. *
  213.       MINWRK = 1
  214.       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  215.          MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
  216.          MINWRK = MAX( 1, 2*N )
  217.          IF( .NOT.WANTVS ) THEN
  218.             MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'SN', N, 1, N, -1 ), 2 )
  219.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'SN', N, 1,
  220.      $          N, -1 ) ) )
  221.             HSWORK = MAX( K*( K+2 ), 2*N )
  222.             MAXWRK = MAX( MAXWRK, HSWORK, 1 )
  223.          ELSE
  224.             MAXWRK = MAX( MAXWRK, N+( N-1 )*
  225.      $               ILAENV( 1, 'ZUNGHR', ' ', N, 1, N, -1 ) )
  226.             MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'SV', N, 1, N, -1 ), 2 )
  227.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'SV', N, 1,
  228.      $          N, -1 ) ) )
  229.             HSWORK = MAX( K*( K+2 ), 2*N )
  230.             MAXWRK = MAX( MAXWRK, HSWORK, 1 )
  231.          END IF
  232.          WORK( 1 ) = MAXWRK
  233.       END IF
  234.       IF( LWORK.LT.MINWRK ) THEN
  235.          INFO = -15
  236.       END IF
  237.       IF( INFO.NE.0 ) THEN
  238.          CALL XERBLA( 'ZGEESX', -INFO )
  239.          RETURN
  240.       END IF
  241. *
  242. *     Quick return if possible
  243. *
  244.       IF( N.EQ.0 ) THEN
  245.          SDIM = 0
  246.          RETURN
  247.       END IF
  248. *
  249. *     Get machine constants
  250. *
  251.       EPS = DLAMCH( 'P' )
  252.       SMLNUM = DLAMCH( 'S' )
  253.       BIGNUM = ONE / SMLNUM
  254.       CALL DLABAD( SMLNUM, BIGNUM )
  255.       SMLNUM = SQRT( SMLNUM ) / EPS
  256.       BIGNUM = ONE / SMLNUM
  257. *
  258. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  259. *
  260.       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
  261.       SCALEA = .FALSE.
  262.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  263.          SCALEA = .TRUE.
  264.          CSCALE = SMLNUM
  265.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  266.          SCALEA = .TRUE.
  267.          CSCALE = BIGNUM
  268.       END IF
  269.       IF( SCALEA )
  270.      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
  271. *
  272. *
  273. *     Permute the matrix to make it more nearly triangular
  274. *     (CWorkspace: none)
  275. *     (RWorkspace: need N)
  276. *
  277.       IBAL = 1
  278.       CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
  279. *
  280. *     Reduce to upper Hessenberg form
  281. *     (CWorkspace: need 2*N, prefer N+N*NB)
  282. *     (RWorkspace: none)
  283. *
  284.       ITAU = 1
  285.       IWRK = N + ITAU
  286.       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
  287.      $             LWORK-IWRK+1, IERR )
  288. *
  289.       IF( WANTVS ) THEN
  290. *
  291. *        Copy Householder vectors to VS
  292. *
  293.          CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
  294. *
  295. *        Generate unitary matrix in VS
  296. *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
  297. *        (RWorkspace: none)
  298. *
  299.          CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
  300.      $                LWORK-IWRK+1, IERR )
  301.       END IF
  302. *
  303.       SDIM = 0
  304. *
  305. *     Perform QR iteration, accumulating Schur vectors in VS if desired
  306. *     (CWorkspace: need 1, prefer HSWORK (see comments) )
  307. *     (RWorkspace: none)
  308. *
  309.       IWRK = ITAU
  310.       CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
  311.      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
  312.       IF( IEVAL.GT.0 )
  313.      $   INFO = IEVAL
  314. *
  315. *     Sort eigenvalues if desired
  316. *
  317.       IF( WANTST .AND. INFO.EQ.0 ) THEN
  318.          IF( SCALEA )
  319.      $      CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
  320.          DO 10 I = 1, N
  321.             BWORK( I ) = SELECT( W( I ) )
  322.    10    CONTINUE
  323. *
  324. *        Reorder eigenvalues, transform Schur vectors, and compute
  325. *        reciprocal condition numbers
  326. *        (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
  327. *                     otherwise, need none )
  328. *        (RWorkspace: none)
  329. *
  330.          CALL ZTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
  331.      $                RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
  332.      $                ICOND )
  333.          IF( .NOT.WANTSN )
  334.      $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
  335.          IF( ICOND.EQ.-14 ) THEN
  336. *
  337. *           Not enough complex workspace
  338. *
  339.             INFO = -15
  340.          END IF
  341.       END IF
  342. *
  343.       IF( WANTVS ) THEN
  344. *
  345. *        Undo balancing
  346. *        (CWorkspace: none)
  347. *        (RWorkspace: need N)
  348. *
  349.          CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
  350.      $                IERR )
  351.       END IF
  352. *
  353.       IF( SCALEA ) THEN
  354. *
  355. *        Undo scaling for the Schur form of A
  356. *
  357.          CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
  358.          CALL ZCOPY( N, A, LDA+1, W, 1 )
  359.          IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
  360.             DUM( 1 ) = RCONDV
  361.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
  362.             RCONDV = DUM( 1 )
  363.          END IF
  364.       END IF
  365. *
  366.       WORK( 1 ) = MAXWRK
  367.       RETURN
  368. *
  369. *     End of ZGEESX
  370. *
  371.       END
  372.